# delimit ; 

************************************ FE*****************************;
**** CHANGE MARCH 2018 - Added 2 year rise @ death for medex graphs and redefined final member deaths as also single households
clear all; 
set mem 300m;	
set more 1 ;  
drop _all;
program drop _all;
capture log close;
set trace on;
global log_medex=1;

cd C:\Dropbox\hrs\wealthcouples;

local data_select=0; * total medex;
*local data_select=1; * medicaid medex;
*local data_select=2; * OOP medex;


if `data_select'==0{;
global graphname medex;
global graphtitle OOP plus Medicaid Expenses;
};
else if `data_select'==1{;
global graphname medicaid;
global graphtitle Medicaid Expenses;
};
else{;
global graphname OOP;
global graphtitle OOP Medical Expenses;
};
log using medexFE2.log, replace ; 


*----------------programs here----------------------------;
*---------------------------------------------------------;
program define fixresid;
			reg res b22-b40 ag11-ag65
ab11-ab66;

predict predres;

fixcoh;
summ predres;
predict meanres; summ predres meanres res;
end;
program define fixcoh;
	version 3.1;
		local i = 22;
		while `i' <=40{;
		replace b`i'=0;
		local i=`i'+1 };

replace b40=1;

fixabag;
end;

*this program generates a polynomial in adist;
program define genPIshift_ln;
        version 3.1;
	replace PI=Pip1;
	makedist;
	gen lmedcostp1=ageshiftvarc+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace lmedcostp1=lmedcostp1+(vage+ PI*bvPI+ PI2*bvPI2  +PIage*bvPIage )/2;
	replace PI=Pip2;
	makedist;
	gen lmedcostp2=ageshiftvarc+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace lmedcostp2=lmedcostp2+(vage+ PI*bvPI+ PI2*bvPI2  +PIage*bvPIage )/2;
	replace PI=Pip3;
	makedist;
	gen lmedcostp3=ageshiftvarc+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace lmedcostp3=lmedcostp3+(vage+ PI*bvPI+ PI2*bvPI2  +PIage*bvPIage )/2;
	replace PI=Pip4;
	makedist;
	gen lmedcostp4=ageshiftvarc+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace lmedcostp4=lmedcostp4+(vage+ PI*bvPI+ PI2*bvPI2  +PIage*bvPIage )/2;

	gen _20th_percentile=exp(lmedcostp1);
	gen _40th_percentile=exp(lmedcostp2);
	gen _60th_percentile=exp(lmedcostp3);
	gen _80th_percentile=exp(lmedcostp4);
drop lmedcostp*;
end;

program define genPIshift_lev;
        version 3.1;
	replace PI=Pip1;
	makedist;
	gen medcostp1=ageshift+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace PI=Pip2;
	makedist;
	gen medcostp2=ageshift+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace PI=Pip3;
	makedist;
	gen medcostp3=ageshift+PI*bPI+PI2*bPI2+(PIage*bPIage);
	replace PI=Pip4;
	makedist;
	gen medcostp4=ageshift+PI*bPI+PI2*bPI2+(PIage*bPIage);
	
	
	gen _20th_percentile=medcostp1;
	gen _40th_percentile=medcostp2;
	gen _60th_percentile=medcostp3;
	gen _80th_percentile=medcostp4;

drop medcostp*;
end;

*this program generates a polynomial in adist;
program define makedist;
        version 3.1;
	drop PI2 PI3 PIage;
	gen PI2=PI*PI;
	gen PI3= PI2*PI;
	gen PIage=PI*age1;
end;



program define genPIshift2_ln;
        version 3.1;

		gen low_couple=b_low+ageshiftvarc +womanhealshift + manhealshift ;
		gen low_male=b_low+ageshiftvarc+single_mshift+vsingle_m/2+manhealshift ;
		gen low_female=b_low+ageshiftvarc+single_wshift+vsingle_w/2 +womanhealshift ;
		
		
		gen high_couple=b_high+ageshiftvarc +womanhealshift+ manhealshift ;
		gen high_male=b_high+ageshiftvarc+single_mshift+vsingle_m/2+manhealshift ;
		gen high_female=b_high+ageshiftvarc+single_wshift+vsingle_w/2 +womanhealshift ;
	
	
	gen _20th_percentile_couple=exp(low_couple);
	gen _20th_percentile_male=exp(low_male);
	gen _20th_percentile_female=exp(low_female);

	gen _80th_percentile_couple=exp(high_couple);
	gen _80th_percentile_male=exp(high_male);
	gen _80th_percentile_female=exp(high_female);

end;

program define genPIshift2_lev;
        version 3.1;

		gen low_couple=b_low+ageshift+womanhealshift +manhealshift ;
		gen low_male=b_low+ageshift+single_mshift +manhealshift ;
		gen low_female=b_low+ageshift+single_wshift +womanhealshift ;
		
		
		gen high_couple=b_high+ageshift+womanhealshift +manhealshift ;
		gen high_male=b_high+ageshift+single_mshift+manhealshift ;
		gen high_female=b_high+ageshift+single_wshift+womanhealshift; 
		
	
	gen _20th_percentile_couple=low_couple;
	gen _20th_percentile_male=low_male;
	gen _20th_percentile_female=low_female;

	gen _80th_percentile_couple=high_couple;
	gen _80th_percentile_male=high_male;
	gen _80th_percentile_female=high_female;

end;


cap program drop drawHS_combo;
program define drawHS_combo;
*This program draws and saves the predicted profiles for a given health status;


if ${log_medex} == 1 {;
*We have already replaced the conditional mean to include the variance in the log specification;
*replace ageshiftvarc=ageshiftvarc+(vage+ PI*bvPI+ PI2*bvPI2 +PIage*bvPIage )/2;

genPIshift2_ln;
};
if ${log_medex} == 0 {;
genPIshift2_lev;
};

if ${log_medex} == 0 {;
gen _20th_percentile_to_male=_20th_percentile_couple;
replace _20th_percentile_to_male=_20th_percentile_male if age>79;
replace _20th_percentile_to_male=_20th_percentile_to_male+bwomandied +low*bwomandiedPI +widower_shift if age==80; 
replace _20th_percentile_to_male=_20th_percentile_to_male+bwomandied +low*bwomandiedPI +widower_shift if age==81; 

gen _20th_percentile_to_female=_20th_percentile_couple;
replace _20th_percentile_to_female=_20th_percentile_female if age>79;
replace _20th_percentile_to_female=_20th_percentile_to_female+bmandied +low*bmandiedPI+widow_shift if age==80; 
replace _20th_percentile_to_female=_20th_percentile_to_female+bmandied +low*bmandiedPI+widow_shift  if age==81;

gen _80th_percentile_to_male=_80th_percentile_couple;
replace _80th_percentile_to_male=_80th_percentile_male if age>79;
replace _80th_percentile_to_male=_80th_percentile_to_male+bwomandied +high*bwomandiedPI+widower_shift if age==80;
replace _80th_percentile_to_male=_80th_percentile_to_male+bwomandied +high*bwomandiedPI+widower_shift if age==81;

gen _80th_percentile_to_female=_80th_percentile_couple;
replace _80th_percentile_to_female=_80th_percentile_female if age>79;
replace _80th_percentile_to_female=_80th_percentile_to_female+bmandied +high*bmandiedPI+widow_shift  if age==80; 
replace _80th_percentile_to_female=_80th_percentile_to_female+bmandied +high*bmandiedPI+widow_shift  if age==81; 

};

if ${log_medex} == 1 {;
gen _20th_percentile_to_male=_20th_percentile_couple;
replace _20th_percentile_to_male=_20th_percentile_male if age>79;
replace _20th_percentile_to_male=_20th_percentile_to_male*exp(bwomandied +low*bwomandiedPI+bvwomandied/2+widow_shift ) if age==80; 
replace _20th_percentile_to_male=_20th_percentile_to_male*exp(bwomandied +low*bwomandiedPI+bvwomandied/2+widow_shift ) if age==81; 
gen _20th_percentile_to_female=_20th_percentile_couple;
replace _20th_percentile_to_female=_20th_percentile_female if age>79;
replace _20th_percentile_to_female=_20th_percentile_to_female*exp(bmandied +low*bmandiedPI+bvmandied/2 +widower_shift ) if age==80; 
replace _20th_percentile_to_female=_20th_percentile_to_female*exp(bmandied +low*bmandiedPI+bvmandied/2 +widower_shift ) if age==81; 
gen _80th_percentile_to_male=_80th_percentile_couple;
replace _80th_percentile_to_male=_80th_percentile_male if age>79;
replace _80th_percentile_to_male=_80th_percentile_to_male*exp(bwomandied +high*bwomandiedPI+bvwomandied/2+widower_shift ) if age==80;
replace _80th_percentile_to_male=_80th_percentile_to_male*exp(bwomandied +high*bwomandiedPI+bvwomandied/2+widower_shift ) if age==81; 
gen _80th_percentile_to_female=_80th_percentile_couple;
replace _80th_percentile_to_female=_80th_percentile_female if age>79;
replace _80th_percentile_to_female=_80th_percentile_to_female*exp(bmandied+high*bmandiedPI+bvmandied/2 +widow_shift ) if age==80; 
replace _80th_percentile_to_female=_80th_percentile_to_female*exp(bmandied+high*bmandiedPI+bvmandied/2 +widow_shift ) if age==81; 
};

replace _80th_percentile_couple=_80th_percentile_couple/1000;
replace _80th_percentile_to_male=_80th_percentile_to_male/1000;
replace _80th_percentile_to_female=_80th_percentile_to_female/1000;

replace _20th_percentile_couple=_20th_percentile_couple/1000;
replace _20th_percentile_to_male=_20th_percentile_to_male/1000;
replace _20th_percentile_to_female=_20th_percentile_to_female/1000;



preserve;
collapse _20th_percentile_to_female _20th_percentile_to_male _20th_percentile_couple _80th_percentile_to_female _80th_percentile_to_male  _80th_percentile_couple, by(age);


la var _20th_percentile_to_female "20th percentile to female";
la var _80th_percentile_to_female "80th percentile to female";
la var _20th_percentile_to_male "20th percentile to male";
la var _80th_percentile_to_male "80th percentile to male"; 
la var _20th_percentile_couple "20th percentile couple";
la var _80th_percentile_couple "80th percentile couple";

local yaxis_opt ylabel(0(5)30);
 local yaxis_opt `yaxis_opt' ymtick(0(2.5)30);

*if;
if "${healstat}"=="NH" { ;

if ${log_medex} == 0 {;
local yaxis_opt ylabel(0(20)100);
local yaxis_opt `yaxis_opt'  ymtick(0(10)100);
};
if ${log_medex} == 1 {;

local yaxis_opt ylabel(0(50)350);
local yaxis_opt `yaxis_opt'  ymtick(0(25)350);
};
};

******************** the key graph;
twoway connected _20th_percentile_to_female _20th_percentile_to_male _20th_percentile_couple age, 
 `yaxis_opt' msymbol(O Oh O) mcolor(red blue green) lcolor(red blue green) lpattern(solid solid solid) 
ti("${graphtitle}", size(medsmall)) sub(" 20th Percentile of Permanent Income, Initially Couples",size(med)) l1(1000s 2014 dollars) ytitle("") saving(medex, replace) legend(order(3 1 2));
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${healstat}_${graphname}_initial_couples20.eps, replace;
graph export graphs/${healstat}_${graphname}_initial_couples20.emf, replace;
*graph export c:\research\hrs\wealthcouples\graphs\medex_initial_couples20.eps, replace;
*graph export c:\research\hrs\wealthcouples\graphs\medex_initial_couples20.emf, replace;

twoway connected _80th_percentile_to_female _80th_percentile_to_male  _80th_percentile_couple age, 
`yaxis_opt' msymbol(S Sh S) mcolor(red blue green) lcolor(red blue green) lpattern(solid solid solid) 
ti("${graphtitle}", size(medsmall)) sub(" 80th Percentile of Permanent Income, Initially Couples",size(med)) l1(1000s 2014 dollars) ytitle("") saving(medex, replace) legend(order(3 1 2));
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${healstat}_${graphname}_initial_couples80.eps, replace;
graph export graphs/${healstat}_${graphname}_initial_couples80.emf, replace;
*graph export c:\research\hrs\wealthcouples\graphs\medex_initial_couples80.eps, replace;
*graph export c:\research\hrs\wealthcouples\graphs\medex_initial_couples80.emf, replace;
sort age;
di "health status: `healstat'";

by age: sum _*;
restore;
drop low_couple low_male low_female high_couple high_male high_female _20th_percentile_* _80th_percentile_*;
end;


set trace off ;


set trace off ;

*------------------------------------------------------------;


use dataprep2_4states;
*use E:\Dropbox\hrs\wealthcouples\dataprep2_4states;
drop *shift*;
drop bPI*;




*group cohort 1 with cohort 2;
replace cohort2=1 if cohort1==1;
replace cohort1=0;

gen cohort8=0;
replace cohort8=1 if cohort==.;
replace cohort=8 if cohort8==1;

* correct concept is household out of pocket + Medicaid;
rename hhmedcost_b medcost_b;
rename hhmedicaid_pay medicaid_pay;
*replace medcost_b=hhmedcost_b;
*replace medicaid_pay=hhmedicaid_pay;
*replace medcost=medcost_b+medicaid_pay;

*select the dependent variable;
if `data_select'==0{; *total medcost;
gen medcost=medcost_b+medicaid_pay;
};
else if `data_select'==1{; *medicaid;
gen medcost=medicaid_pay;
};
else{; *OOP medcost;
gen medcost=medcost_b;
};

*replace medcost_b=1000 if medcost_b<1000 & medcost_b~=.;

sum medcost, det;
bys manheal womanheal: sum medcost, det;

egen tempmean=mean(medcost), by(manheal womanheal);
replace medcost=0.2*tempmean if medcost<0.2*tempmean & medcost~=.;
drop tempmean;

tab manheal womanheal,m;

*replace manheal=0 if hhstatus==3;
*replace womanheal=0 if hhstatus==2;
*tab manheal womanheal,m;

gen manbadheal=(manheal==2);
gen mannursheal=(manheal==1);
*gen mandiedheal=(manheal==0);
gen womanbadheal=(womanheal==2);
gen womannursheal=(womanheal==1);
*gen womandiedheal=(womanheal==0);
*replace medcost=. if mandied==1|womandied==1;

sort HHID wave;
gen lmanbadheal=manbadheal[_n-1] if HHID==HHID[_n-1] & wave <wave[_n-1] + 3;
gen lmannursheal=mannursheal[_n-1] if HHID==HHID[_n-1] & wave <wave[_n-1] + 3;

gen lwomanbadheal=womanbadheal[_n-1] if HHID==HHID[_n-1] & wave <wave[_n-1] + 3;
gen lwomannursheal=womannursheal[_n-1] if HHID==HHID[_n-1] & wave <wave[_n-1] + 3;

* new medcost variable here;
sum medcost;

*No bottom coding of the sum here, instead we bottom code hhmedcost_b in dataprep2;

* log medcost;
gen lmedcost=ln(medcost);

sum lmedcost;

gen tempE=exp(r(mean)+r(var)/2);
sum tempE;


reg lmedcost;
gen meanmed=_b[_cons];
gen residmed=lmedcost-meanmed;
gen residsq=residmed*residmed;
reg residsq;
gen varpred=_b[_cons];
gen predmean=exp(meanmed+varpred/2);
gen pred975=exp(meanmed+2*sqrt(varpred));
sum predmean pred975 medcost, d;

* BREAK



*drop if lmedcost==.;
*replace heal=1 if dead==1; * in dataprep2 i do this if died==1.  changing it to "if dead==1" gets at medical expenses of those who died in previous waves;
*drop if heal==.; * drops 3 obs;
*sum *medc* PI male heal age1;

/*
drop age age1 age2 age3 age4 age5;
gen age=hhage;
gen age1=hhage;
gen age2=age1^2;
gen age3=age1^3;
gen age4=age1^4;
gen age5=age1^5;
*/

*ALTEERNATIVE SPECIFICATION IN THIS TEST FILE;
replace single_m=1 if hhstatus==0 & mandied==1 & hhstatus[_n-1]~=3; *& hhstatus[_n-1]==2;
replace single_w=1 if hhstatus==0 & womandied==1 & hhstatus[_n-1]~=3; *& hhstatus[_n-1]==1;

sort age1;
by age1: sum medcost lmedcost;

*gen manhealage=manheal*age1;
*gen womanhealage=womanheal*age1;

drop couple;
gen couple=0;
replace couple=1 if hhstat==1;

gen coupleage=couple*age1;
gen coupleage2=couple*age2;

gen single_mage=single_m*age1;
gen single_wage=single_w*age1;

gen single_mage2=single_m*age2;
gen single_wage2=single_w*age2;

gen mandiedPI=mandied*PI;
gen womandiedPI=womandied*PI;


*construct NH interactions;
gen single_mnursheal=single_m*mannursheal;
gen single_wnursheal=single_w*womannursheal;

gen single_mlnursheal=single_m*lmannursheal;
gen single_wlnursheal=single_w*lwomannursheal;


*tab medcost;

sort wave;
by wave: sum age1 medcost;

makedist;




* OLS VERSUS FIXED EFFECTS;
* OLS is not consistent with what JOHN, Cristina and I have decided on--it is only descriptive;
if FE==0{;
reg lmedcost age1 age2 age3  heal healage maleage PIage;

reg lmedcost age1 age2 age3  heal healage maleage PIage PI PI2 male;
gen  ageshiftOLS=_b[_cons]+(_b[age1]*age1)+(_b[age2]*age2) +(_b[age3]*age3);
sort age1;
by age1: sum ageshift*;
};

if FE==1{;
* fixed effects ;
*xtreg lmedcost age1 age2 age3 age4 age5, fe i(HHID);
*xtreg lmedcost age1 age2 age3 age4 age5 heal, fe i(HHID);
if ${log_medex} == 1 {;
sum lmedcost age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal single_mage single_wage single_mage2 single_wage2 mandied womandied mandiedPI womandiedPI PIage single_wnursheal single_mnursheal;

xtreg lmedcost age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal single_mage single_wage single_mage2 single_wage2 mandied womandied mandiedPI womandiedPI PIage single_wnursheal single_mnursheal single_wlnursheal single_mlnursheal , fe i(HHID);
*xtreg lmedcost age1 single_m single_w manheal womanheal  single_mage single_wage PIage , fe i(HHID);
*manhealage womanhealage;
predict lmedcostp;


};
if ${log_medex} == 0 {;
xtreg medcost age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal  single_mage single_wage single_mage2 single_wage2 mandied womandied mandiedPI womandiedPI PIage single_wnursheal single_mnursheal single_wlnursheal single_mlnursheal, fe i(HHID);
*xtreg lmedcost age1 single_m single_w manheal womanheal  single_mage single_wage PIage , fe i(HHID);
*manhealage womanhealage;
predict medcostp;
};

* get predictions;
gen ageshift=_b[_cons]+(_b[age1]*age1)+(_b[age2]*age2) +(_b[age3]*age3) +(_b[age4]*age4) ;
gen manbadhealshift=_b[manbadheal] +_b[lmanbadheal]; //+_b[manhealage]*age1; *RORY HAS CHANGED;
gen womanbadhealshift=_b[womanbadheal]+_b[lwomanbadheal]; //+_b[womanhealage]*age1;  *RORY HAS CHANGED;

gen mannurshealshift=_b[mannursheal]+_b[lmannursheal]; //+_b[manhealage]*age1; *RORY HAS CHANGED;
gen womannurshealshift=_b[womannursheal]+_b[lwomannursheal];

gen smannurshealshift=_b[single_mnursheal]+_b[single_mlnursheal];
gen swomannurshealshift=_b[single_wnursheal]+_b[single_wlnursheal];

gen single_mshift=_b[single_m]+_b[single_mage]*age1+_b[single_mage2]*age2;
gen single_wshift=_b[single_w]+_b[single_wage]*age1+_b[single_wage2]*age2;


gen PIshift=_b[PIage]*age1;
gen bPIage=_b[PIage];

gen bmandied=_b[mandied];
gen bwomandied=_b[womandied];

gen bmandiedPI=_b[mandiedPI];
gen bwomandiedPI=_b[womandiedPI];
local bsavelist b bwomandiedPI bmandiedPI bwomandied bmandied bPIage bPI bPI2;
*Added to pick up the constant: RM;
gen b=_b[_cons];

foreach varname in age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal single_mage single_wage single_mage2 single_wage2  single_wnursheal single_mnursheal single_wlnursheal single_mlnursheal{;
gen b`varname'=_b[`varname'];
local bsavelist `bsavelist' b`varname';
};
*close the foreach;
sum b*;

if ${log_medex} == 1 {;
drop lmedcostp; *RORY HAS CHANGED below;
gen lmedcostp=ageshift+(blmanbadheal*lmanbadheal)+(blwomanbadheal*lwomanbadheal)+ (blmannursheal*lmannursheal)+(blwomannursheal*lwomannursheal)+(bmanbadheal*manbadheal)+(bwomanbadheal*womanbadheal)+ (bmannursheal*mannursheal)+(bwomannursheal*womannursheal)+ (single_mshift*single_m)+(single_wshift*single_w)+(PIshift*PI)+mandied*(bmandied+bmandiedPI*PI)+womandied*(bwomandied+bwomandiedPI*PI) + single_wnursheal*bsingle_wnursheal +single_mnursheal*bsingle_mnursheal +  single_wlnursheal*bsingle_wlnursheal +single_mlnursheal*bsingle_mlnursheal;
gen lmedcostres=lmedcost-lmedcostp;
};
if ${log_medex} == 0 {;
drop medcostp; *RORY HAS CHANGED below;
gen medcostp=ageshift+(blmanbadheal*lmanbadheal)+(blwomanbadheal*lwomanbadheal)+ (blmannursheal*lmannursheal)+(blwomannursheal*lwomannursheal)+(bmanbadheal*manbadheal)+(bwomanbadheal*womanbadheal)+ (bmannursheal*mannursheal)+(bwomannursheal*womannursheal)+(single_mshift*single_m)+(single_wshift*single_w)+(PIshift*PI)+mandied*(bmandied+bmandiedPI*PI)+womandied*(bwomandied+bwomandiedPI*PI) +  single_wnursheal*bsingle_wnursheal +single_mnursheal*bsingle_mnursheal +  single_wlnursheal*bsingle_wlnursheal +single_mlnursheal*bsingle_mlnursheal;
gen medcostres=medcost-medcostp;
};


*sum lmedcost*;
* now run an OLS regression on the residuals;
if ${log_medex} == 1 {;
reg lmedcostres PI PI2 cohort5 cohort3 cohort4 cohort6 cohort7 cohort8;
};
if ${log_medex} == 0 {;
reg medcostres PI PI2 cohort5 cohort3 cohort4 cohort6 cohort7 cohort8;
};
predict respred;
* set cohort dummies to 0;
replace cohort2=0;
replace cohort3=0;
replace cohort4=0;
replace cohort6=0;
replace cohort7=0;
replace cohort8=0; //need to regress again after?;
predict meanres;

if ${log_medex} == 1 {;
gen lmedcadj=lmedcost+meanres-respred;
gen medcadj=exp(lmedcadj);
};

save medresid, replace;
*save C:\research\hrs\wealthcouples\medresid, replace;

*Added to pick up the constant: RM;

replace b=b+_b[_cons];
replace ageshift=ageshift+_b[_cons];
gen ageshiftFINAL=ageshift;
gen bPI =_b[PI];
gen bPI2=_b[PI2];

replace PIshift=PIshift+_b[PI]; * do not add in quadratic term, as we keep track of that separately;
* see if cohort effects are important to worry about (as it turns out, they are);


* figure out variance of residual;
if ${log_medex} == 1 {;
gen resres=lmedcostres-respred;
};
if ${log_medex} == 0 {;
gen resres=medcostres-respred;
};
*drop if resres==.;
egen sdres=sd(resres);
gen varc=(sdres*sdres);

* figure out the unconditional variance;
egen sdmedc=sd(lmedcost);
gen varcmed=sdmedc*sdmedc;

* R2 measure;
gen R2meas=1-varc/varcmed;

*gen varc2=lmedcostres*lmedcostres; //do we need this?;
gen varc1=resres*resres;
reg varc1;
reg varc1 age1 ;
reg varc1 age1 age2 PI;



*xtreg varc1 age1 age2 age3 heal healage PIage, fe i(HHID);
*gen vageFE=_b[_cons]+(_b[age1]*age1)+(_b[age2]*age2) +(_b[age3]*age3);




reg varc1 age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal PI PI2 PIage mandied womandied single_wnursheal single_mnursheal single_wlnursheal single_mlnursheal; //single_mage single_wage ;
*reg varc1 age1 single_m single_w manheal womanheal PI PI2 PIage ; //single_mage single_wage ;
*reg varc1 age1 age2 age3 age4 single_m single_w  PI PI2 PIage cohort2 cohort3 cohort4  cohort6 cohort7 cohort8;

predict varc1p;
//gen vage=_b[_cons]+(_b[age1]*age1);//+(_b[age2]*age2) +(_b[age3]*age3) +(_b[age4]*age4);
gen vage=_b[_cons]+(_b[age1]*age1)+(_b[age2]*age2) +(_b[age3]*age3) +(_b[age4]*age4);

gen vsingle_m=_b[single_m];//+(_b[single_mage]*age1);
gen vsingle_w=_b[single_w];//+(_b[single_wage]*age1);

gen vmanbadheal=_b[manbadheal]+ _b[lmanbadheal];
gen vmannursheal=_b[mannursheal] + _b[lmannursheal];
gen vwomanbadheal=_b[womanbadheal] + _b[lwomanbadheal];
gen vwomannursheal=_b[womannursheal] + _b[lwomannursheal];

gen vsmannursheal=_b[single_mnursheal] + _b[single_mlnursheal];
gen vswomannursheal=_b[single_wnursheal] + _b[single_wlnursheal];

gen vPI=_b[PI]+(_b[PIage]*age1);
gen vPI2=_b[PI2];

gen bv=_b[_cons];
gen bvPI=_b[PI];
gen bvPIage=_b[PIage];
gen bvPI2=_b[PI2];
local bvsavelist bv bvPI bvPIage bvPI2;

foreach varname in age1 age2 age3 age4 single_m single_w lmanbadheal lwomanbadheal lmannursheal lwomannursheal manbadheal womanbadheal mannursheal womannursheal mandied womandied single_wnursheal single_mnursheal single_wlnursheal single_mlnursheal{;
gen bv`varname'=_b[`varname'];
local bvsavelist `bvsavelist' bv`varname';
};
*close loop;

* sort age1;
* by age1: sum vage*;
*two conn vage vageFE age1;

sum varc* R2meas;



gen ageshiftvarc=ageshift;

* gen medcostp=exp(lmedcostp); * this term is meaningless, in the aggregate;
};


*rename age1 age;
sort age HHID;

*for the next 10 lines we generate the normalised residuals to estimate the AR1 and transitory decomposition;
preserve;

gen normalizedresid=resres/sqrt(varc1p);
*demean by wave - no assymptotic difference (asympt mean=0);
sum wave;
local maximumwave=r(max);

forvalues waveloop=3(1)`maximumwave'{;

sum normalizedresid if wave==`waveloop', meanonly;

local temp_mean=r(mean);
replace  normalizedresid=normalizedresid-`temp_mean' if wave==`waveloop';
};

*keep the variables needed;

keep normalizedresid wave HHID;

save normalizedresid, replace;

restore;
*Sort by health status and sum predictions vs data:;
if ${log_medex} == 0 {;
*resres is the gap between the prediction (including PI effects) and the observed data;
gen exponentres=resres;
sum exponentres, det;
bys manheal womanheal: sum exponentres, det;
};
if ${log_medex} == 1 {;
*expectation of medcost:;
gen tempexponent=exp(lmedcostp+respred+varc1p/2);
*total error;
gen exponentres=medcost-tempexponent;
gen res95=medcost-exp(lmedcostp+respred+1.65*sqrt(varc1p));
gen res50=medcost-exp(lmedcostp+respred);
sum exponentres, det;
bys manheal womanheal: sum exponentres res95 res50, det;

sum exponentres res95 res50, det;
sum exponentres res95 res50 if manheal==1|womanheal==1, det;
};



*BREAK;



drop if age>100|age<70;
drop if age==age[_n-1];
save ${graphname}prof, replace;
*save c:\research\hrs\wealthcouples\medexprof, replace;
*drop if age>100;

***********************************;
* GRAPHS                           ;
***********************************;

set scheme s1mono;



* COUPLE, BOTH GOOD HEALTH;

*In order to take off a blue or pink background delete "plotregion(ifcolor(pink*0.08))" in a line below;
if ${log_medex} == 1 {;
genPIshift_ln;

};
if ${log_medex} == 0 {;
genPIshift_lev;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);
twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Couple in Good Health", size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}gc.eps, replace;
graph export graphs/${graphname}gc.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgc.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgc.emf, replace;
**********numbers behind the graphs****************;
restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age ageshiftvarc vage;
drop _*;

* COUPLE, BOTH BAD HEALTH;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+manbadhealshift+womanbadhealshift+vmanbadheal/2+vwomanbadheal/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-manbadhealshift-womanbadhealshift-vmanbadheal/2-vwomanbadheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+manbadhealshift+womanbadhealshift;
genPIshift_lev;
replace ageshift =ageshift-manbadhealshift-womanbadhealshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Couple in Bad Health",size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}bc.eps, replace;
graph export graphs/${graphname}bc.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbc.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbc.emf, replace;
**********numbers behind the graphs****************; restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age ageshiftvarc vage;
drop _*;


* COUPLE, BOTH NURSING HOME;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+mannurshealshift+womannurshealshift+vmannursheal/2+vwomannursheal/2;
genPIshift_ln;
;
replace ageshiftvarc =ageshiftvarc-mannurshealshift-womannurshealshift-vmannursheal/2-vwomannursheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+mannurshealshift+womannurshealshift;
genPIshift_lev;
replace ageshift =ageshift-mannurshealshift-womannurshealshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(40000)200000) ymtick(0(20000)200000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Couple in Nursing Home",size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}nc.eps, replace;
graph export graphs/${graphname}nc.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnc.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnc.emf, replace;
**********numbers behind the graphs****************; restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age ageshiftvarc vage;
drop _*;


* WOMEN, GOOD HEALTH;

*In order to take off a blue or pink background delete "plotregion(ifcolor(pink*0.08))" in a line below;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_wshift+vsingle_w/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-single_wshift-vsingle_w/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_wshift;
genPIshift_lev;
replace ageshift =ageshift-single_wshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid longdash shortdash longdash_dot) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Women in Good Health", size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}gw.eps, replace;
graph export graphs/${graphname}gw.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgw.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgw.emf, replace;
**********numbers behind the graphs****************; restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;


* WOMEN, BAD HEALTH;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_wshift+womanbadhealshift+vsingle_w/2+vwomanbadheal/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-single_wshift-womanbadhealshift-vsingle_w/2-vwomanbadheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_wshift+womanbadhealshift;
genPIshift_lev;
replace ageshift =ageshift-single_wshift-womanbadhealshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Women in Bad Health",size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}bw.eps, replace;
graph export graphs/${graphname}bw.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbw.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbw.emf, replace;
**********numbers behind the graphs****************; restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;


* WOMEN, NURSING HOME;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_wshift+womannurshealshift+vsingle_w/2+vwomannursheal/2 +swomannurshealshift +vswomannursheal/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-single_wshift-womannurshealshift-vsingle_w/2-vwomannursheal/2 -swomannurshealshift -vswomannursheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_wshift+womannurshealshift +swomannurshealshift;
genPIshift_lev;
replace ageshift =ageshift-single_wshift-womannurshealshift-swomannurshealshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(40000)200000) ymtick(0(20000)200000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Women in Nursing Home",size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}nw.eps, replace;
graph export graphs/${graphname}nw.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnw.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnw.emf, replace;
**********numbers behind the graphs****************;
restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;


* MEN, GOOD HEALTH;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_mshift+vsingle_m/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-single_mshift-vsingle_m/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_mshift;
genPIshift_lev;
replace ageshift =ageshift-single_mshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Men in Good Health", size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}gm.eps, replace;
graph export graphs/${graphname}gm.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgm.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexgm.emf, replace;
**********numbers behind the graphs****************;
restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;


* MEN, BAD HEALTH;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_mshift+manbadhealshift+vsingle_m/2+vmanbadheal/2;
genPIshift_ln;

replace ageshiftvarc =ageshiftvarc-single_mshift-manbadhealshift-vsingle_m/2-vmanbadheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_mshift+manbadhealshift;
genPIshift_lev;
replace ageshift =ageshift-single_mshift-manbadhealshift;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(5000)30000) ymtick(0(2500)30000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Men in Bad Health", size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}bm.eps, replace;
graph export graphs/${graphname}bm.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbm.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexbm.emf, replace;
**********numbers behind the graphs****************;
restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;

* MEN, NURSING HOME;
if ${log_medex} == 1 {;
replace ageshiftvarc =ageshiftvarc+single_mshift+mannurshealshift+vsingle_m/2+vmannursheal/2 +smannurshealshift +vsmannursheal/2;
genPIshift_ln;
replace ageshiftvarc =ageshiftvarc-single_mshift-mannurshealshift-vsingle_m/2-vmannursheal/2 -smannurshealshift -vsmannursheal/2;
};
if ${log_medex} == 0 {;
replace ageshift =ageshift+single_mshift+mannurshealshift+smannurshealshift ;
genPIshift_lev;
replace ageshift =ageshift-single_mshift-mannurshealshift-smannurshealshift ;
};
preserve;
collapse _20th_percentile _40th_percentile _60th_percentile _80th_percentile, by(age);

twoway connected _20th_percentile _40th_percentile _60th_percentile _80th_percentile age, legend(off) ylabel(0(40000)200000) ymtick(0(20000)200000) msymbol(T S D O) lpattern(solid solid solid solid) ti("${graphtitle}", size(medsmall)) sub(" by Permanent Income Percentile, Men in Nursing Home", size(medsmall)) l1(2014 dollars) ytitle("") saving(medex, replace);
graph use medex.gph;
graph display, ysize(6) xsize(8.25);
graph export graphs/${graphname}nm.eps, replace;
graph export graphs/${graphname}nm.emf, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnm.eps, replace;
*graph export C:\research\hrs\wealthcouples\graphs\medexnm.emf, replace;
**********numbers behind the graphs****************;
restore;
sort age;
by age: sum _20th_percentile _40th_percentile _60th_percentile _80th_percentile age;
drop _*;


gen low=.2;
gen low2=low^2;
gen high=.8;
gen high2=high^2;
gen lowage=low*age1;
gen highage=high*age1;

*calculate the conditional mean from regression results;
gen b_low=bPI*low+bPI2*low2+(lowage*bPIage);
gen b_high=bPI*high+bPI2*high2+(highage*bPIage);


if ${log_medex} == 1 {;
replace b_low=b_low+(vage+ low*bvPI+ low2*bvPI2 +lowage*bvPIage )/2;
replace b_high=b_high+(vage+ high*bvPI+ high2*bvPI2 +highage*bvPIage )/2;

*select the health status to produce the graph for:;
*GOOD HEALTH!;
*define globals used by PIshift programs:;
gen manhealshift=0;
gen womanhealshift=0;
*When the man dies assign the lagged health status:;
gen widow_shift=0;
*When the woman dies assign the lagged health status:;
gen widower_shift=0;

global healstat GH;

drawHS_combo;
*BH case;
replace manhealshift=bmanbadheal+blmanbadheal+(bvmanbadheal+bvlmanbadheal)/2;
replace womanhealshift=bwomanbadheal+blwomanbadheal+(bvwomanbadheal+bvlwomanbadheal)/2;
*When the man dies assign the lagged health status:;
replace widow_shift=blmanbadheal+(bvlmanbadheal)/2;

*When the woman dies assign the lagged health status:;
replace widower_shift=blwomanbadheal+(bvlwomanbadheal)/2;
global healstat BH;

drawHS_combo;



*NH case;
replace manhealshift=bmannursheal+blmannursheal+(bvmannursheal+bvlmannursheal)/2 ;
replace womanhealshift=bwomannursheal+blwomannursheal+(bvwomannursheal+bvlwomannursheal)/2;
*When the man dies assign the lagged health status - and remmeber that they are now current single so need the single*lag spouse health;
replace widow_shift=blmannursheal+(bvlmannursheal)/2 +bsingle_mlnursheal +bvsingle_mlnursheal/2;

replace single_mshift=smannurshealshift +single_mshift +vsmannursheal/2;;
replace single_wshift=swomannurshealshift +single_wshift +vswomannursheal/2;;

*When the woman dies assign the lagged health status - and remmeber that they are now current single so need the single*lag spouse health:;
replace widower_shift=blwomannursheal+(bvlwomannursheal)/2 +bsingle_wlnursheal +bvsingle_wlnursheal/2;
global healstat NH;
drawHS_combo;
};
if ${log_medex} == 0 {;

*select the health status to produce the graph for:;
*GOOD HEALTH!;
*define globals used by PIshift programs:;
gen manhealshift=0;
gen womanhealshift=0;
*When the man dies assign the lagged health status:;
gen widow_shift=0;
*When the woman dies assign the lagged health status:;
gen widower_shift=0;

global healstat GH;

drawHS_combo;

*BAD HEALTH;
*define globals used by PIshift programs:;
replace manhealshift=bmanbadheal+blmanbadheal;
replace womanhealshift=bwomanbadheal+blwomanbadheal;
*When the man dies assign the lagged health status:;
replace widow_shift=blmanbadheal;

*When the woman dies assign the lagged health status:;
replace widower_shift=blwomanbadheal;
global healstat BH;
drawHS_combo;


*NURSING HOME;
*define globals used by PIshift programs:;
replace manhealshift=bmannursheal+blmannursheal;
replace womanhealshift=bwomannursheal+blwomannursheal;
*When the man dies assign the lagged health status:  - and remmeber that they are now current single so need the single*lag spouse health;
replace widow_shift=blmannursheal+bsingle_mlnursheal;

*When the woman dies assign the lagged health status:  - and remmeber that they are now current single so need the single*lag spouse health;
replace widower_shift=blwomannursheal+bsingle_wlnursheal;
replace single_mshift=smannurshealshift+single_mshift;
replace single_wshift=swomannurshealshift +single_wshift;
global healstat NH;
drawHS_combo;

};
*;


replace single_mshift=-bsingle_mnursheal+single_mshift;
replace single_wshift=-bsingle_wnursheal+single_wshift;

drop widow_shift widower_shift manhealshift womanhealshift;
*Sort by health status and sum predictions vs data:;


use ${graphname}prof, replace;
*use c:\research\hrs\wealthcouples\medexprof, replace;

keep `bsavelist' `bvsavelist';
/*keep age ageshiftFINAL manbadhealshift  mannurshealshift single_mshift womanbadhealshift womannurshealshift single_wshift PIshift bPI2 
vage  vmanbadheal vmannursheal vwomanbadheal vwomannursheal  vsingle_m  vsingle_w 
 vPI vPI2;
order age ageshiftFINAL manbadhealshift  mannurshealshift single_mshift womanbadhealshift womannurshealshift single_wshift PIshift bPI2 
vage  vmanbadheal vmannursheal vwomanbadheal vwomannursheal  vsingle_m  vsingle_w 
 vPI vPI2;
*/
order `bsavelist' `bvsavelist';

keep if _n==1;
sum;

keep 
b
bage1
bage2
bage3
bage4
bPI
bPI2
bPIage
bsingle_m
bsingle_w
bmanbadheal
bwomanbadheal
bmannursheal
bwomannursheal
blmanbadheal
blwomanbadheal
blmannursheal
blwomannursheal
bsingle_mage
bsingle_wage
bsingle_mage2
bsingle_wage2
bmandied
bwomandied
bmandiedPI
bwomandiedPI
bv
bvage1
bvage2
bvage3
bvage4
bvPI
bvPI2
bvPIage
bvsingle_m
bvsingle_w
bvmanbadheal
bvwomanbadheal
bvmannursheal
bvwomannursheal
bvmandied
bvwomandied
bvlmanbadheal 
bvlwomanbadheal 
bvlmannursheal 
bvlwomannursheal 
bsingle_wnursheal 
bsingle_mnursheal
bsingle_wlnursheal
bsingle_mlnursheal
bvsingle_wnursheal 
bvsingle_mnursheal
bvsingle_wlnursheal
bvsingle_mlnursheal;


order
b
bage1
bage2
bage3
bage4
bPI
bPI2
bPIage
bsingle_m
bsingle_w
bmanbadheal
bwomanbadheal
bmannursheal
bwomannursheal
blmanbadheal
blwomanbadheal
blmannursheal
blwomannursheal
bsingle_mage
bsingle_wage
bsingle_mage2
bsingle_wage2
bmandied
bwomandied
bmandiedPI
bwomandiedPI
bsingle_mnursheal
bsingle_wnursheal
bsingle_mlnursheal
bsingle_wlnursheal
bv
bvage1
bvage2
bvage3
bvage4
bvPI
bvPI2
bvPIage
bvsingle_m
bvsingle_w
bvmanbadheal
bvwomanbadheal
bvmannursheal
bvwomannursheal
bvmandied
bvwomandied
bvlmanbadheal 
bvlwomanbadheal 
bvlmannursheal 
bvlwomannursheal 
bvsingle_mnursheal
bvsingle_wnursheal
bvsingle_mlnursheal
bvsingle_wlnursheal;

export excel using ${graphname}prof_FE, sheetreplace firstrow(variables);


outsheet using ${graphname}prof_FE,nonames replace;
*outsheet using c:\research\hrs\wealthcouples\medexprof_FE,nonames replace;
drop _all;
program drop _all;
log close;
